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Abstract 

A new method for the computation of electromagnetic scattering from arbitrary two- 
dimensional bodies is presented. The method combines the finite element and boundary 
element methods leading to a system for solution via the Conjugate Gradient FFT_ al- 
gorithm. Two forms of boundaries aimed at reducing the storage requirement of the 
boundary integral are investigated. It is shown that the boundary integral becomes con- 
volutional when a circular enclosure is chosen, resulting in reduced storage requirement 
when the system is solved via the Conjugate Gradient FFT method. The same holds for 
the ogival enclosure, except that some of the boundary integrals are not convolutional 
and must be carefully treated to maintain the 0(N ) memory requirement. Results for 
several circular and ogival structures are presented and shown to be in excellent agree- 
ment with those obtained by traditional methods. 
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Chapter 1 

Introduction 


Many methods exist for the numerical solution of two-dimensional (2-D) scattering 
problems. Moment methods have traditionally dominated the frequency domain solu- 
tion approaches though more recently, hybrid finite element methods have gained much 
popularity. The relative simplicity in treating complex arbitrary composite structures 
is a major reason for this. Also, the resulting system of equations in a finite element 
implementation is sparse and banded leading to a low ( 0(N )) storage required for large 
scale applications. 

To formulate the hybrid finite element method for scattering computations, the struc- 
ture is enclosed in a ficticious boundary. Within the boundary, the finite element method 
is used to solve a weak representation of the Helmholtz equation and further, to satisfy 
the radiation condition, an approximate absorbing boundary condition (ABC) [1] may 
be placed on the ficticious boundary. The ABC's axe popular because they result in 
a banded sub-matrix structure. However, they require additional unknowns since the 
enclosure must be placed at a distance approximating the far field region. An alternative 
to the ABCs is to match the fields within the enclosure to an eigenfunction expansion 
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(unimoment method) [2] or to employ the boundary integral equation [3, 4, 6]. The 
unimoment method produces a dense square sub-matrix with dimension proportional to 
the number of modes. It also requires the truncation of an infinite series which may be 
slowly convergent for irregular structures thus resulting in a large storage requirement. 

Previously, the authors introduced a method [6] which resulted in an 0(N ) storage 
requirement. By choosing a rectangular ficticious boundary, some of the integrals in 
the boundary integral equation become convolutions amenable to evaluation via the fast 
Fourier transform (FFT). Provided the conjugate gradient (CG) algorithm is used for 
the solution of the system, the discretization and evaluation of the convolution integrals 
requires only 0(N) storage (while the remaining “cross-terms” must be stored in an 
efficient manner). Another important feature is that the order of the FFT need only 
be applied on the ficticious boundary making it preferable to the traditional CGFFT 
method, which requires the evaluation of FFTs that have the same dimensionality as 
that of the structure. 

It is possible to choose other boundaries that result in convolutional integrals, and in 
this report we consider circular and ogival enclosures. Clearly, a circular enclosure would 
be attractive for circular scatterers whereas an ogival boundary will be more attractive 
for those structures conforming to this boundary. In the case of the circular boundary 
the entire integral is convolutional ensuring the 0(N ) memory demand of the system 
provided an iterative solution is used. When an ogival enclosure is used the integral 
becomes convolutional only if the observation and source points are on the same arc, but 

an efficient storage scheme is again required for the remaining “cross-terms” l . 

1 “cross terms” refer to integrals for which the source and observation points are not on the same arc 
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The hybrid finite element methods presented herein will be referred to as the com- 
bined finite element- boundary element methods (BE/FE-CGFFT), because the source 
and observation points associated with the boundary integral shaxe a common contour 
as is the case with the traditional boundary element method [8]. In contrast, for the 
combined finite element - boundary integral formulation (FE-CGFFT) discussed in [6] 
the source and observation points are on adjacent contours. This choice of a single or 
double contour is, of course, arbitrary and only for convenience since in either case, the 
same number of unknowns are required for either case. 

In the following sections, the pertinent BE/FE-CGFFT formulations are developed 
for the circular and ogival boundaries. Results for several circular and ogival structures 
are presented and shown to be in excellent agreement with that obtained by traditional 
methods. 
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Chapter 2 

Analysis 


Consider the plane wave 1 

$ inc (p) = z4 nc (p) = ze^ cos (*-*°> (2.1) 

illuminating a composite cylinder as shown in Fig. 2.1 and we axe interested in computing 
the scattered field. For the application of the Finite Element - Boundary Element Method 
the target is enclosed in a ficticious circular or ogival boundary as shown in Figs. 2.2 and 
2.3. Within the boundary the finite element method is used to solve the Helmholtz 
equation 


V • [v(p)V<t>(p)} + klv(p)4>{p) = 0 (2.2) 

where 

4>(p) = E z (p), u{p) = — v(p) = € r (p) (2.3) 

PryP) 

for E-polarization and 

4>{p) = H,(p), u(p) = v(p) = p r (p) (2.4) 

5 An e ,Mt time convention is assumed and suppressed. 
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Figure 2.1: Geometry of the scatterer 
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for H-polarization. The free-space wave number is k Q — Uy/fi 0 e 0 and fi r and e r are the 
relative permeablility and permittivity, respectively. On the boundary T a the Helmholtz 
integral equation 


m = <P' nC (p) - / r {G{p,-p a ) [£- a <t>(Pa)] - * 0 .) [-£- a G (P'P«i\ } dl « ( 2 ‘ 5 ) 

provides the required boundary constraint, implicitly satisfing the radiation condition. 
In (2.5) 

G(p,p a ) = -^H^ 2) (k 0 \p - p a |) (2.6) 

is the 2-D free space Green’s function where Ho 2 \>) denotes the zeroth order Hankel 
function of the second kind. Also, denotes differentiation with respect to the outward 
normal, whereas ~p and p 0 are the usual source and observation points, respectively, and 

IP “Pal = \/(z-Za) 2 + (2/-&>) 2 (2.7) 
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Figure 2.2: Partially discretized body in a circular enclosure 

2.1 Case 1: Circular Enclosure 

2.1.1 Discretization of the Scatterer and Field Quantities 

The region enclosed by T a , denoted as i2 a , is discretized into N c finite elements as 
illustrated in Fig. 2.2. In the figure, p a is the radius of the circle and a a is the integration 
angle along this boundary (Further definitions for the finite element mesh are indicated 
in Table 2.1, while the definitions of the field vectors are indicated in Table 2.2.). We 
note that nodes along T a are equispaced with angular displacement A. 

2.1.2 Derivation of the Finite Element Matrix 

The weighted residual expression over each element may be written [9] 

J I Rci'N'dQ' = 0 i = 1,2,3 (2.8) 

n e 
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Definitions for Finite Element Mesh 

N n — total number of nodes in the finite element mesh 

N g = total number of unknowns 

N e = total number of elements in the finite element mesh 
jV 0 = number of nodes or elements on T a 
Nd = number of nodes or elements on 

Table 2.1: Definitions for the finite element mesh 


Definitions of Field Vectors (in terms of field unknowns at nodal points) 

4> a = fields at the nodes on T a 

4 > a = normal field derivatives at the nodes on T* 

<f)j = fields corresponding to region I enclosed by T 0 and T 4 
<f>d = fields at the nodes on the conductor boundary 

Table 2.2: Definitions of the field vectors 
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where 


R U‘ = - J- 


u(x,y)j^d> e (x,y) 


d_ 

dy 


3 

u(x,y)—4> e ( x ,y)\ - k 2 0 v(x,y)<p e {x,y) 


(2.9) 


and 


4>~4> = Y.P 

e=l 


( 2 . 10 ) 


where <f> e is zero outside element e . In (2.9) is an approximation to 4> the eth 
element, and N f is the ith shape function associated with the eth element. Substituting 
(2.9) into (2.8) and invoking the divergence theorem yields 

+ j d9 e 

[ N^ e dT e = Q (2.11) 

Jr c 



—u 
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d dN? d4> e dNf 


dx dx dy dy 


where T e denotes the contour enclosing the eth element. Additionally, 


d<t> e 


* - 


( 2 . 12 ) 


is zero outside element e. Summing over N e elements we obtain 

+ 


E 



dtf_ dNf df_ dNf 
dx dx dy dy 


+ 


iv4>'N'\ dQ.' 

N* . N d 

E / WPdTl + E / . N&tti = 0 

3 = 1 Jr * 3=1 JV d 


(2.13) 


where the summations over s refer to the elements with sides adjacent to the ficticious 
(T a ) and conducting (I\j) boundaries. The integral over the conducting boundary van- 
ishes and if no conductor is present, then Tj is not present. When <f> = H z , the normal 
derivative of the field is zero on the conductor and the field unknowns on the boundary 
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are allowed to “float” (i.e., the boundary condition is “naturally” satisfied). Finally, 
when <j> = E z , imposing the Dirichlet condition during assembly of the finite element 
system results in the elimination of those equations associated with the integral over IV 
Proceeding with the discretization, the field and its derivative within each element 
may be expanded into a linear combination of shape functions 

& = i:m ( 2 - 14 ) 

:=i 

r = E wit (2.i5) 

*:=i 

Substituting (2.14) and (2.15) into (2.11) we obtain 

3 3 

= ° ( 2 - 16 ) 

j=l Ar=l 

where 



dN]_dN]_ dNJ_ dNf 

dx dx dy dy 


- k*vN?N, 


:} 


d(l e 


(2.17) 


and 


b’k = [ N’ k N*dr* a 

J r» 


For linear triangular elements, Nf are given by 


N? = W '(a? + b'x + cly) 


with 




1 , 

- det 

2 


1 y\ 
i A 

1 X% l tl 



(2.18) 


(2.19) 


(2.20) 
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( 2 . 21 ) 


_ e „ e _ .e _e .,e 

a i — x j Vk X kVj 

b e i = (2.22) 

Ci = xl-x) (2.23) 


and (xf , 2 /f ) being the coordinates of the ith node of the eth element. From (2.19) 


dNf 

dx 

II 

(2.24) 

dNf 

dy 

ii 

(2.25) 


Using (2.24), (2.25) and the identity 


a 


plq\ 


*("0W^ = 2«' ()1 + f + 2)! 


(2.26) 


aij in (2.17) reduces to 


?/ c O e 

/ _ ( k<* ue i „c^e\ _ 


«<; = 4 ^m + ^-ky-d + s^) 


(2.27) 


where 




1 if i — j 


(2.28) 


0 otherwise 

We note that in deriving (2.27) we have assumed that u and v (the reciprocal of the 
material constitutive parameters) are constant within each element and are given by u e 
and t? e , respectively. 

To find an algebraic expression for b 9 lk , we may reparametize the integral in (2.18) as 


b 


3 

ik 


r 

J a i 


ai +A 


P’ P£r a da 


(2.29) 


11 



where Pf and P? are given by 


Integrating, we have 




! _ Q ~ Q i 
a 2 “ a l 

a — a{ 
a 2 - <*’ 


% = r -f-(i<k + 1 ) 


(2.30) 

(2.31) 


(2.32) 


Substituting the previous equations into (2.13) a sparse matrix is obtained for the 
nodal fields that has the form 


A a a A a j 0 —B a a 

: 


(j>a 


0 

Ma All Aid 0 


4>i 


0 

0 Adi Add 0 


4>d 


0 

0 0 0 0 


'Pa 


0 


(2.33) 


In this, the values of the elements in the submatrix A pq are the contributions associated 
with the nodes in group (region or boundary) p which are connected directly to the nodes 
in group q. Also, 


N a 


[•5ao]»'fc — ^ Kk ~ + 46, fc + 6,+i fc) 


5 = 1 


(2.34) 


The last row in (2.33) has been intentionally left blank to imply a need for another set 
of equations relating the fields and its derivatives on T a . This additional set of equations 
is produced by discretizing the boundary integral equation. 
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The last row in (2.33) has been intentionally left blank to imply a need for another set 
of equations relating the fields and its derivatives on T a . This additional set of equations 
is produced by discretizing the boundary integral equation. 


2.1.3 Evaluation of the Boundary Integral 

The boundary integral in (2.5) may be rewritten in cylindrical coordinates via the 
transformations 

|/3 — P a | = \x(pcosa - p a cosa a ) - y(psina - /> a sina a )| 

= \[p 2 + ~pl - 2 ~pp a cos(a - a a ) (2.35) 

where {p y ot) and (p a ,a a ) are the usual source and observation points in cylindrical co- 


ordinates. For \p\ = |p a |, 

\p~Pa\ = 2p|sin(a=a*)| (2.36) 

and the Green’s function and its normal derivative may be written as 

G{p,? a ) = -taw (2k oPa sin (^)) (2.37) 

J~G(p,? a ) = J -^H\ 2) (2k oPa sin(^))sm(^). (2.38) 

on a 4 

We may now write (2.5) as [8] 

^<f>(p, a) = <t>' nc (p, a) - fo(p, a) + fi(p, a) (2.39) 

where as a result of (2.37) and (2.38) 

fo{p,cx) = -{pa / i>{p a ,a a )Hi 2) (2kop a sm(z=2*-))da a (2.40) 

4 Jo 
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(2.41) 


fl(p, a ) = ^ Pa f Q ^(PayOta)!!^ (2k 0 p a sin(^^-))sm(^^)da a 


with 


'J’(Pa,a a ) = —4>(p a ,a a ) 
O pa 


(2.42) 


The factor of ~ in (2.39) accounts for the singularity associated with #^(*) and the -f 
in (2.40) and (2.41) denotes principle value. 

We may now discretize (2.40) by expanding the field using pulse basis functions as 




TpiPa) ^ ^ 
i=i 


(2.43) 


P&(<*a ~ Otj) = { 


(2.44) 


where 

1 \ct a - a_j| < A 
0 otherwise 

and A is the angular width of the integration cell as indicated in Fig. 2,2. Thus, the 
discrete version of (2.40) may be written as 

fo(p,a) = f 3 + lrt 2) WoPaSm(2=^))da a (2.45) 

q j = i J «>- T 

Performing point collocation and letting u' = a — a a , we have 

^ f (a,-a 7 H4 


fo{p,oti) = -f J 3 #i 2) (2k 0 p a sin (£)) cfa' 

4 ; = 1 J 

which may be written in compact form as 

fo(p,t*i) = - ] ~y £ - Of j ) 

where 

/■(“.-<*>)+$• ... 

fto(a, -a J )= f Hy> {2k 0 p a sin (f))du 

J (a, -a, )- ■* 


(2.46) 


(2.47) 


j=l 


(or.-aj)-T 


(2.48) 
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It is clear that (2.47) is in the form of a discrete convolution and can thus be written as 


fo(p,a) = DFT- 1 { DFT( 0) • DFT(h 0 )} 


(2.49) 


where the elements of ho are given by 

- i] } p = 0 

r(P+£) A 


h 0 (pA) = { 


(2.50) 


l f(p- })A H ° ] { 2k oPaSin (t)) du ' p = N « - 1 

and 7 ~ 1.781 [10]. Through a similar analysis, the field may be approximated by the 
expansion 


Na 

^(Pa> a a) — ^ ^ ^ > a( < ^o ““ )^j 

J = i 

and by substutiting this into (2.41), we obtain 

fi {Pi <*• ) = - «i) 

;=i 

where 

/ix(ai - a j) = -f ‘ 2 (2fc 0 p a sin (^)) sin (*£) <fu' 

d (a.-orjJ-A 

Clearly, (2.52) may again be written in operator form as 

/i(p,a) = DFT~ l {DFT(j>) . DFTihr)} 


where 


fti(pA) = 


Wa -sinf) p=0 

/(?JL)A ^1 (2) ( 2fc oPa sin sin (^) du' p = 1, ..., N a - 


Point matching (2.39) at each node results in the system 

1& = C C - /o(p, a.) + /i(p,«i) 


(2.51) 


(2.52) 


(2.53) 


(2.54) 


(2.55) 


(2.56) 
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which may be written in operator form as 

M aa 4>a - Ua^a = 4%* 

where 

[i-l« = 

A final system is obtained by combining (2.57) with (2.33) to yield 


0 B aa 


■ 

i 


0 

Aia An Aid 0 


4>i 


0 

0 Adi Add 0 


4>d 


0 

M aa 0 0 ~~L aa 


, 

4 

j 


_ 


which can be solved via the conjugate gradient algorithm to obtain the nodal fields. 


(2.57) 

(2.58) 

(2.59) 
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Definitions for Finite Element Mesh 

N n = total number of nodes in the finite element mesh 

N g = total number of unknowns 

N e = total number of elements in the finite element mesh 
N a = N ai + N a2 = number of nodes on T a 

r a = r ai + r a2 

Table 2.3: Definitions for the finite element mesh 

2.2 Case 2: Ogival Enclosure 

2.2.1 Discretization of the Scatterer and Field Quantities 

The region within T a , denoted ii a , is discretized into N e finite elements and a partial 
discretization is shown in Fig. 2.3 for the circular case. In Fig. 2.3, 

A p = angular displacement between nodes on T 0p 
p ap = radii of 

a Qp = angular integration variable along F ap 
t — distance between centers of curvature of T 0 

Up 

Vc p = y - coordinate of the center of curvature of T 0p 

Further definitions for the finite element mesh are indicated in Table 2.3, and the field 
vector definitions are indicated in Table 2.4. 
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Definitions of Field Vectors (in terms of field unknowns at nodal points) 


<f>a p = fields corresponding to the nodes on T ap , p = 1,2 

V>a p = fields corresponding to the midpoints of the nodes on T ap 

< Pa p = fields at the nodal midpoints on T ap 

<j>l = fields corresponding to region I enclosed by T a and T d 

4> d = fields corresponding to the nodes on the T d 



Table 2.4: Definitions of the field vectors 
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2.2.2 Derivation of the Finite Element Matrix 


The derivation of the finite element matrix follows that described in section 2.1.1 
with the exception of the matrix B aa . Consider the ogival boundary as indicated in fig. 
2.3. The boundary contour T a is comprised of two arcs labeled T ai and r a2 » which form 
the vertices of the ogive where they meet. At the vertices the unknown normal field is 
discontinuous and will therefore be evaluated at the midpoint. Also, in evaluating the 
contour integral, the field derivative will be expanded in terms of pulse basis functions, 
rather than linear functions. This results in a different B aa matrix and involves the 
replacement of Pj in (2.29) by the pulse basis function expansion 


Fa(« ~ otj) = 


1 if 0 < \a-aj\ < f 


0 otherwise 


By integrating in cylindrical coordinates we then obtain 


(2.61) 


bij — 0 ^i,j+ 1)> 3 — I* * — 1*2 


(2.62) 


where l e is the length of the eth boundary element along T a and is equal to /? ap A p for 
Tap? P = 1,2. Performing a summation over all boundary elements then yields 

[^aa]jj = ^ ^ij = + £«,j+ l) (2.63) 

c=l Z 

where l J is the length of the j th element since the j th w node” (associated with the 
unknown ipj) is at the center of the j th boundary element. 

The remainder of finite element analysis for this case proceeds exactly as in section 

2 . 1 . 2 . 
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2.2.3 Evaluation of the Boundary Integral 

The evaluation of the boundary integral along an ogival contour is similar to that 
described for the circular boundary. For integration and observation points on the same 
arc, the integrals become convolutions. On the other hand, when the integration and 
observation points reside on different contours, the integrals have no special form and 
must be discretized and stored in memory as efficiently as possible. 

The distance between the source and observation points in terms of cylindrical coor- 
dinates for points on the same arc is given by 

IP - PaJ = \[P 2 + Pa p - 2ppa p COs(<* - Q ap ) p= 1,2 (2.64) 

When the source and observation points are along different arcs, (2.64) becomes 

Ip, - Po p l = \j (p cos - Pa cos a ap ) 2 + (p sin a, - p a sin a 0j> + y c<J - y^) 2 

P,q= 1,2 (2.65) 


in which the subscript a p refers to the integration coordinates along contour p and the 
subscript q refers to the observation coordinates. Also, y ^ is the y-coordinate of the 
center of curvature for contour p for p = 1,2. For further reference we note that (2.65) 
may be also rewritten a- s 



2 


- Pa 2 I 


1 


> 



- 2pi p a2 cos(ai - a 0j ) + t 2 T 2t(pi sinai - p a2 sina a2 ) (2.66) 

2 1 2 i 2 2 1 1 


where t = - y Ci . 
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To discretize (2.5), the fields are expanded as 


N ai 

4>(pa,Ota) ^ E 
J = 1 


- + Y p a( a “ - otj)<Pj 

j=N ai +1 


N aj 

rp{p a ,a a ) ~ YL P ±( a ° 

j=\ 


N a 

~ <*])$] + Y P A(«a 

j=No, +1 


- oj )h 


where as before 


(^a ^j) — 


1 if |a a - otj\ < & 
0 otherwise 


and 


<Pj = 2 ^; + 4>j+ 1) 


Substutiting (2.67), (2.68) and (2.69) into (2.5) then yields 


1 f - r<*>+A i 

■z4>(Pi,oti) = 4> >nc (Pi,oti) — ^ Y'bi T Go{P\,p ai ,a\ - a ai )p ai da ai 
z L i=i y 


N 0 


„ fOj + A-2 

+ / , t G 0 (pi,p a2 , a i, a a2 )p a2 da a2 

j=N ai +l J a ) 


ro 

-Ew/ 

• t ^ a 


i . / k ° f j +Al <9 

j_l ^ afj ^Pai 


^o(Pl * Pd\ y ^1 & ax') P a\ a\ 


^ rot ]+&2 Q 'y 

^ ^ fij T ^ C”o(Pl i Pa2 > *-*1 » ^*<12 )Pa2 dot^ r 

J= AT a +1 y «J a ^2 J 


when the observation point is on T ai and 


1 f i - yaj+Ai 

-T0(/} 2 ,Q 2 ) = <p inC {P2y<*2) ~ f G 0 {p2,Pa x ^^2yOL ax )p { 

1 K j = l J a > 

Ng r a 

+ E */ 

j=iV ai +l 


ai da ai 


fdj + Ai 


Go(,P2i Pcl2 ) )Pfl2 


aj 


(2.67) 


(2.68) 


(2.69) 


(2.70) 


(2.71) 
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N« 


a f a Q 

~ ^2 &j T (*o{p 2 y Pa\ > a 2> ^ai )pa\ dot ai 


U ^ j *, **, 

ya J+ A 2 Q 

" L WJ 

i=N a , +1 7 "■> 




^-*o(P2j Pa 2 ? ^-2 ^-a 2 )P<* 2 j (2. 


72) 


for observation on r a2 . Performing point collocation at the nodal midpoints, (2,71) and 
(2.72) further reduce to 

1 ~ roti-aj + % 

O0(pi,^i + i) = ^ mc (pi»<*i + i) - \ zZ'tPj f G 0 (pi,p ai ,u)p ai du 

1 2 2 ^ j= i </«*-«;+$■ 


N a 


f C*j -fA 2 


' Ta ^ f c 

+ ^2 M G 0 (^l,^>a 2? «t+^ a O2)Pa 2 ^a 2 

j=N ai +1 ^ “* 

^ rai-a 7 +£ $ 

— ^2 <fj T A ;J 7 {?o(pi*P<li ? W )pOi du 

j= 1 J «•-»>-■ 


J+£ ^P«l 


^ ^ /ar>+A 2 £ 'l 

— } " T ■q” ^*o(/* 1 * /*a 2 5 a i + l. ? tt a 2 )Pa 2 da a2 ^ (2.73) 

i=AT 0l +l 7 ap ° 2 J 


for observation on r oi and 


1 f * - /’ a J" h ^ 1 

«^(P2i«i + i) = ^ nC (P2,a t -+l)- G 0 (p 2 ,Pa 1 ,0( |+ l,^a 1 )p ) 

^ 2 7 ^ j = i J *) 7 

^ * Qi - Qj + A 

+ 12 ^ T Go(p 2 ,Pa 2 ,u)p { 

jsAT.,+1 ^ °*-«> + T 


ai ^*ai 


du 


Na 


. raj+Ai ^ 

— /2 T ;7” Go(p2i Pa\ y » **ai )Poi do: aj 

i=X y «> 


^Pai 

Na />ar l -a J +£ P 'I 

- £ <Pj f .j— G 0 (p 2 ,p aj ,u)p a2 du\ (2.74) 

J= ^+i J *,-<,,+$ dp* 7 J 

where the in the subscript refers to the ficticious “node" midway between the actual 
nodes. A system of equations can now be obtained by testing (2.73) and (2.74) at a 
sequence of points on the contours. This yields 


1 


-{ 


1 


7^1031 = <ra" C ~ { ^110ai + A20a 2 “ TjA/jiCx 4 > ( 




- 2^12^2^02 ) 
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-C 2 <fid 2 — y>'a 2 C — 1^21^0! + ^22^a 3 ~ ^ ~ 2^ 22 ^ 2< ^ a 2 | (2.75) 


which can be alternatively written as 


M\\C\ Q12C2 

Q21C1 M.22C2 


D 


Ln Pi 2 
P21 L22 


(fiai 

< t > a 2 

0 ai 

^a 2 


vs 


V 


a 2 


(2.76) 


In this C is a matrix with I’s down the main and super diagonals, and <p are the value 
of the incident field evaluated at the nodal midpoints. The matrix D accounts for the 
double use of the nodes at the endpoints and the remaining elements are given by 


Qpq = ~~^Qpq 


(2.77) 

(2.78) 

(2.79) 


for p = 1,2 in which 


and 


[Afpp]- ■ 


= / 
= / 


£*,-£*,+ 




')T J 


Go^Pp t Pa p > U)p ap du 


0,-0; + - 


0,-0; + ^- $Pa 

/•o.+A, 


Goi^Ppi Pap » ri)papd\L 


[Lpp\,j 

fa, +fl, 

= / G 0 (Pp, Pa q ,0‘i + L, ( Xaq)Pa q dO( aq 

Jat, 1 


/■Oj+A, £ 

= ~Q~p G 0 {Pp)Pa q > a ,+l > a a q )Pa q dot aq 


(2.80) 

(2.81) 

(2.82) 


(2.83) 
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More explicitly, upon evaluation of the integrals 


[Mpp]- 




= < 


i^pplij — 


\ Ql i = h£\k£T±) 

1 IlJij 4 J aj (fc oV /- V <") 


= < 


j^T*- [ko pa p {^f ~ sin i = j 

J a '~ a ’+£ H[ 2) (2k oPap sin f )sin fd« i ? j 

ai-aj + -f- r it *! 

i = j 


f a ' "'^4 Ho 1] {2k 0 p^s\n'4)du i^j 

or. - otj + -j* 


(2.84) 


(2.85) 


p a2 — pi cos(ai — a aj ) ± t sin a, 


“2 


1 21 J »j 4 i 


da a2 (2.86) 


(2.87) 


where the upper sign corresponds to the upper set of subscripts on P or Q, while the 
lower sign corresponts to the lower set of subscripts. Introducing the definitions 


K, = 


K 2 = 


D 


MuCi Q12C 2 

Q21C1 AI22C2 

•^n A2 
P21 7/22 

the system (2.76) may be combined with that derived via the finite element method to 
obtain 


( 2 . 88 ) 


(2.89) 


4aa 4 a / 0 B aa 


Ala All Aid 0 

0 Adi Add 0 

A', 0 0 K 2 


" 

" 


- 


4*a 


0 


<t>[ 


0 


<Pd 


0 


1 

1 




(2.90) 
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We note that (2.90) can be solved via the CG algorithm to take advantage of the con- 
volution operators M and L in reducing the memory requirements. This algorithm is 
given next. 
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2.3 A CGFFT Algorithm 


The CG algorithm to be employed for solving a system of the form 



II 

(2.91) 

as in (2.60) and (2.90) is as 

follows: 


Initialize the residual and search vectors 


ib 

= II [ <f> inc 0 0 0 ] r 111=11 b 111 

(2.92) 

s 

= 

(2.93) 

r (l) 

= b — s 


s 

= A a r (1) 

(2.94) 

Is 

= II * 111 

(2.95) 

0(0) 

= 77 1 

(2.96) 

pW 

= 

(2.97) 

Iterate for k = 


3 = ApW 

(2.98) 


IM <N 
®0 

II 

•> 

(2.99) 


« W = 77 1 

(2.100) 


*<*+»> = ^) + a W p W 

(2.101) 


r (*+i) = r (k) _ a (k) p (k) 

(2.102) 


7r = || r (.k+l) ||1 

(2.103) 
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3 


= A a r^ +1 ) 


Is 

,(*+i) 


*113 


Is 1 


P'-.-, - p (k) + p(k) s (k) 

Terminate when k = N g or < tolerance. 


{^} = {-SB/} + {sfe}- 


where 


and 


{SBI} = 


( 5 fe ) = 


Ai a An Aid 0 

0 0 Add 0 

0 3 0 0 

For the adjoint operations, we have 

0 0 0 


0 0 0 0 

0 0 0 0 

0 0 0 0 

Ki 0 0 K 2 

;fl A a I 0 J9q 


T 

" 


Z\ 


Z* 


z 3 


ZA 


T 

" 


z\ 


Z 2 


Z3 


Z 4 


{*b/} = 


0 0 0 0 
0 0 0 0 


0 0 0 iq 


*1 

" 


Z\ 


Z2 


z 3 


z 4 


(2.104) 

(2.105) 

(2.106) 

(2.107) 

(2.108) 


(2.109) 


( 2 . 110 ) 


( 2 . 111 ) 
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and 


{sfe} = 


. 4 a 

aa 

Ki 

0 

0 


Z\ 

A a 

A Ia 

Vi 

Ah 

0 


z 2 

0 

Aa 

A dl 

Ah 

0 


%3 

Baa 

0 

0 

0 


Z4 


( 2 . 112 ) 


In each case, the operation is performed such that only the resulting vector {5} needs 


to be stored. 
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Chapter 3 

Scattered Field Computation 


In this section the expressions for the scattered field and radar cross-section are 
developed for both the circular and ogival boundaries. The scattered fields may be 
computed from the identity 


= - / { g ( 7 , a ) [ 3 ^®.)] - } Jr - P-i) 


and the echowidth is then found from 



3.1 Circular Boundary 


The scattered field expression (3.1) may be written as 


0 s (p,a) = -f 0 (p,a) + fi(p, a) 


(3.3) 


where 


fo(p,a) = -^p a J o xl){pa,a a )H^ ) (fco^V + Pa - 2pp a cos (a - a a )) da a (3.4) 


and 


flip, a) = 


j ft* H x (2) (i k 0 s/p 2 + Pa - 2 pp a cos (a - a a )J 

Tpa*o / <^(pa,a 0 ) / * == 7 ? (p a -pcos(a - a a ))da a { 3.5) 

4 7o sjp 2 + pi - 2pp a cos(a - a a ) 

To evaluate the integrals in (3.4) and (3.5) we invoke the field expansions (2.43) and 
(2.51). We have 

* f a J + ^ 


/o(p,Q) = -tp a J 1 ^ H { 0 2) (fco^/p 2 + Pa - 2 /V»a cos ( Q - “a)) rfa a (3-6) 


and 


/i (/>.“) = 


j £• . /■«, + $• ff} 2) (wp 2 + Pa - 2ppl cos(a - o«)) 

\/p 2 + P 2 - 2pp a cos(a — a„) 


(p a - pcos(a - a a )) dad(3.7) 


where the remaining integrals over the subsections must be evaluated numerically for ar- 
bitrary observation. However, for far-held computations (p — » oo), the Hankel functions 


may be approximated as 


HV\kp) ~ 



-i*p 


(3.8) 
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and since 


\Jp 2 + P\~ 2 pp a cos(a - a a ) ~ 
(3.6) and (3.7) become 


p for amplitude 

p — p a cos(a — a a ) for phase terms 


and 


4 V wk 0 p j9l 
Substutiting (3.10) and (3.11) into (3.3) we obtain 




(3.9) 


/o(p,a) = (3.10) 


flip , a) = \Hf- e ~ ik0P E h cos ( a ~ a i )e^“ COB <"-^) (3.11) 


N a Na 

+ * 0 £^ cos ( a - 

. ;=i j=i 


4 V TT^oP 

and from (3.2) the echowidth is given by yields the echowidth 

0>qA) 2 


3.12) 


A ax 8?r 


N a N a 

j he jkop ° + k 0 2 h cos (a - aj )e jkopa 

>=i i=i 


(3.13) 
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3.2 Ogival Boundary 


Following the same discretization scheme used in Section 2.2.3, (3.1) may be written 
as 

, N »l ^ iVa 

<t> 3 (p,a) = Y'ipf u ( P ,a,a J ) + V V>/i 2 (p, a, aj) 

^ ; = 1 j=M>, +i 

N “i _ /V a 

- 53 - 53 4>f22(p,ot,a } ) > (3.14) 

J = 1 i=^Vaj +i J 

where 


flp(p1 a 1 a ]) — 

f a :+b p 

I G 0 (p, p ap , a, Q! ap )Pa p doi ap 

Ja tj 

(3.15) 

f2p(p,Ot,Otj) = 

ya ; +Ap Q 

1 n (*o(P> Pdp yG"> **a p )/*a p 

«/ar ; t'Pop 

(3.16) 



(3.17) 


in which 


G 0 {p,Pa p ,a,ot ap ) 

- (fro\/p 2 + P2 p - 2 ~pp7 p cos (a - a 0p ) + - 2y Cp (psina - p ap sina 0p ))(3.18) 

and 

d 

G 0 {p,p ap ,a,a ap ) 

__ j'fc 0 + “ 2p/5 °p cos ( Q ~ a “p) + y =p “ 2 yc„(psina - p ap sina ap )j 

4 + Pa p - 2 PPa P cos(a - a ap ) + J/cp “ 2j/ Cp (p sin a - p ap sin a a J 

[Pa p ~ pc os(a - a 0p ) + y Cp sin a ap ](3.19) 
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and are the corresponding ^-coordinates of the arc T ap . Using the large argument 
approximation for the Hankel functions 




(3.20) 


and the approximation 


\Jp + p\~ 2pPa p cos(a - a ap ) + + 22 / Cp (psina - p 0p sina ap ) 


= < 


for amplitude terms 


(3.21) 


P ~ Pa p cos(a - a ap ) + y Cp sin a for phase terms 
for p — > oo, the Hankel function simplifies to 


Ho 2 H\fy^7)) ~ ^J-^-e~ jkoP e~ jk °\.~ Pa p cos (“-“«p)-»c p sin a] 


(3.22) 


Similarly, 
gi (2) (/M) 


^2^2. [. ycv ] ~ -jj ~JL e -jl'°Pe - jk °[- p °’> ™«] cos(a - a ap X3,23) 

\J-Vcj,-) V 

Substutiting these into (3.15) and (3.16) and performing midpoint integration yields 

/ lp (p,a,Oj) = A [-<■•» (3.24) 

4 y 7rA; 0 p 

4 y 


cos(aj + - a 0p ) (3.25) 


Thus, from (3.14) 
<t>ff(P 


) JV 

, a) = -^y^« -ifcoP { - jAip., £ cos(a- 0 2-^) 


~y Cl sin o | 




-jA 2 p aj £ ^->*.1-^ «*(«—,-¥) 

i=N«,+ 1 


-y C2 1 
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original page is 

OF POOR QUALITY 


N ai ^ 

—k 0 A\p ax Y c«(a-«, -%>■)_ fcl sina,] cos(a . + ) 


i=i 


N 

- k 0 A 2Pa 2 Y 4>e-M-^ c ^ a - a x-^)-yc 2 sina] cos (a J + ^-a a2 )} (3.26) 

■ ii * J 


Na 

I 

j=AT 01 +l 

and by substituting this into (3.2) yields the echowidth 

iv t 


A 8tt 


N ai 


1 jAi/Ja, yj xPe~ jko l~ Pa x sina] 

i=l 

y\r 

+jA 2Pa2 Y ^ e -jk 0 [- Pa 2 c<*( a - aj -%-)-y C 2 smc,] 


+ k 0 Ai Pai Y &~ jk ° [ - pa ' c “ ( “-^-^ L )-Vc 1 sina] cos(a . + Al_ _ aoi } 


j=l 


M, 


+ koA 2Pa2 Y co 8 (a-a,-AM-yc 2 sina] cos(q . + ^ ) 

J=iV«, +1 ^ 


(3.27) 
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Chapter 4 

Results 


The scattering patterns of several circular and ogival cylinders for both E- and H- 
polarization are shown in the figures to follow. Figs. 4. 1-4.3 contain circular geome- 
tries both coated and uncoated, while Figs. 4. 4-4. 6 contain coated and uncoated ogival 
structures. The echowidth is computed for each structure and compared to the results 
of the series solution for the circular bodies and moment method [5, 4] for the ogival 
structures. As seen in all cases, the generated patterns via the hybrid BE/FE-CGFFT 
formulation are in excellent agreement with the corresponding data based on the Mei 
Series and Moment Method Solutions. 
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Perfectly Conducting Cylinder R=.5 k (E-pol) 



Perfectly Conducting Cylinder R=.5 k (H-pol) 



Figure 4.1: E z and II Z bistatic echowidth from a perfectly conducting circular cylinder 
of radius 0.5A. 
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Coated Conducting Cylinder R=.5 X (E-pol) 



Angle [deg] 


Coated Conducting Cylinder R=.5 X (H-pol) 



Figure 4.2: E z and H z bistatic echowidth from a perfectly conducting circular cylin- 
der with a conductor radius of .5A and a coating thickness of .05A containing material 
properties c T = 5 - j 5, = 1.5 - jO.5. 
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Coated Conducting Cylinder R=3 A (E-pol) 



Coated Conducting Cylinder R=3 A (H-pol) 



Angle [deg] 


Figure 4.3: E z and Hz bistatic echowidth from a coated circular cylinder with a conductor 
radius of 3A and coating thickness of 0.05A with material properties ( T = 5 — j 5, /r r = 
1.5 -jO. 5. 
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Echo WidthA [dB] Echo WidthA [dB] 


.5 x 1 A PEC Ogive (E-pol) 



Angle [deg] 


.5 x 1 A PEC Ogive (H-pol) 



Figure 4.4: E z and II Z backscatter echowidth from a 0.5 x 1A perfectly conducting ogive. 
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.5 x 1\ Coated Conducting Ogive (E-pol) 



Angle [deg] 


.5 x IX Coated Conducting Ogive (H-pol) 



Figure 4.5: E z and H z backscatter echowidth from a .5 X 1A perfectly conducting ogive 
with a 0.05A thick material coating containing the properties c r = 3-j5,/x r = 1.5- jO. 5. 
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1 x 4 A. Coated Conducting Ogive (E-pol) 



Angle [deg] 


1 x 4A. Coated Conducting Ogive (H-pol) 



Angle [deg] 


Figure 4.6: E z and H z backscatter echowidth from a 1 x 4A perfectly conducting ogive 
with a .05A thick material coating containing the properties e r = 3- j5,fi r = 1.5- j0.5. 
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Chapter 5 

Conclusions and Future Work 


The scattering from targets surrounded by ogival and circular boundaries has been 
presented. The finite element method produces the usual sparse sub-matrix, while a 
discrete version of the boundary integral results in a dense sub-matrix. The mathematical 
boundary enclosing the scattering structure may be judiciously chosen such that the 
boundary integrals are convolutional. As a result, they become amenable to evaluation 
via the FFT and leads to an 0(N ) storage requirement. Among the circular and ogival 
boundaries considered, the circular boundary satisfies the above requirements. The 
ogival boundary results in convolutions only when the source and observation points are 
along the same arc, while the non-convolutional cross-terms must be stored efficiently 
to guarantee the required storage requirement. When considering circular and ogival 
structures, the associated circular and ogival boundaries are usually conformal to the 
structure, thus providing an additional reduction in the number of unknowns. 

To validate the method and associated computer code, scattering patterns of several 
circular and ogival structures were given and compared with data generated by proven 
methods. 
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Figure 5.1: Three-dimensional finite-cylinder enclosure. 

The fundamental purpose of this work was to explore formulations that lead to 0(N) 
storage requirements when employed for three-dimensional simulations. The presented 
hybrid technique may prove useful when the surrounding boundary is chosen such that 
most of the terms of the boundary integral are convolutional. The remaining “cross- 
terms” however, must be stored efficiently and some sort of interpolation should be used 
for their evaluation. Based on this study, a suitable three-dimensional enclosure is a 
truncated circular cylinder as illustrated in Fig. 5.1. 
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